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ABSTRACT 



Context. Knowledge about the coronal magnetic field is important to the understanding of many phenomena, such as flares and 
coronal mass ejections. Routine measurements of the solar magnetic field vector are traditionally carried out in the photosphere. We 
compute the field in the higher layers of the solar atmosphere from the measured photospheric field under the assumption that the 
corona is force-free. However, those measured data are inconsistent with the above force-free assumption. Therefore, one has to apply 
some transformations to these data before nonlinear force-free extrapolation codes can be applied. 

Aims. Extrapolation codes of cartesian geometry for medelling the magnetic field in the corona do not take the curvature of the 
Sun's surface into account. Here we develop a method for nonlinear force-free coronal magnetic field medelling and preprocessing of 
photospheric vector magnetograms in spherical geometry using the optimization procedure. 

Methods. We describe a newly developed code for the extrapolation of nonlinear force-free coronal magnetic fields in spherical 
coordinates over a restricted area of the Sun. The program uses measured vector magnetograms on the solar photosphere as input and 
solves the force-free equations in the solar corona. We develop a preprocessing procedure in spherical geometry to drive the observed 
non-force-free data towards suitable boundary conditions for a force-free extrapolation. 

Results. We test the code with the help of a semi-analytic solution and assess the quality of our reconstruction qualitatively by 
magnetic field line plots and quantitatively with a number of comparison metrics for different boundary conditions. The reconstructed 
fields from the lower boundary data with the weighting function are in good agreement with the original reference fields. We added 
artificial noise to the boundary conditions and tested the code with and without preprocessing. The preprocessing recovered all main 
structures of the magnetogram and removed small-scale noise. The main test was to extrapolate from the noisy photospheric vector 
magnetogram with and without preprocessing. The preprocessing was found to significantly improve the agreement between the 
extrapolated and the exact field. 
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1. Introduction 

The magnetic field in the solar corona dominates over non- 
magnetic forces s uch as plasm a pressure and gravity because of 
low plasma beta (lGaryll200ll) . Knowledge of the coronal mag- 
netic field is therefore important in understanding the structure 
of the coronal plasma and obtaining insights into dynamical 
processes such as flares and coronal mass ejections. Routine 
measurements of the solar magnetic field are still mainly car- 
ried out in the photosphere. Therefore, one has to infer the 
field strength in the higher layers of the solar atmosphere 
from the measured photospheric field based on the assump- 
tion that the corona is force-free. The extrapolation methods 
invol ved in this assumption inclu de potential field extrapola- 
tion dSchmidtl '1964'; 'Semel' '1967), linear force-free field ex- 



ing for several reasons. Mathematically, problems regarding the 
existence and uniqueness of various boundary value problems 
dealing with nonl inear force-free fields remain unsolved (see 
lAmari et"aLll2006l for details). Another issue is their numerical 
analysis of given boundary values. An additional complication 
is to derive the boundary data from observed photospheric vec- 
tor magnetic field measurements, which are consistent with the 
force-free assumption. High noise in the transverse components 
of the measured field vector, ambiguities regarding the field di- 
rection, and non-magnetic forces in the photosphere complicate 
the task of deriving suitable boundary conditions from measured 
data. For a more complete review of existing methods for com- 
puting nonlinear for ce-free coronal m agnetic fi elds, we refer t o 



trapol ation jChiu & Hiltora,1977- Seehafeil [T978l 119821: ISemell 
119881: IClegg et aL l2000h. and nonlinear force-free field extrap- 



the review papers bylAmari et al. (1997), Schri iver et al.l (|2006|) . 
iMetcalfet al.r(l2008 



and lWiegelniand (120081) . 



olation ("Am ari et al.lll997l 1 19991 120061: ICuperman et al.l 1 199 it 



Dem oulin et 
1996; Sakurai 



198lt IValqri etaLl l2005t 



;t all ll^^^L IZUUbt ICupermaii et all 11^^ It 
I992I; iMikic & McClvmontlll994l: iRoumeliotisI 



'Wheatland' 



Wiegelmannl2004l:IWu et al.lll990l:lYan &"Sakural2000.) 



12001 



Among these, the nonlinear force-free field has the most re- 
alistic description of the coronal magnetic field. The computa- 
tion of nonlinear force-free fields is however, more challeng- 



The magnetic field is not force-free in either the photo- 
sphere or the lower chromosphere (with the possible exception 
of sunspot areas, where the field is strongest). Furthermore, mea- 
surement errors, in particular for the transverse field components 
(eg. perpendicular to the line of sight of the observer), would de- 
stroy the compatibility of a magnetogram with the condition of 
being force-free. One way to ease these pro blems is to prepro- 
cess th e magnetograph data as suggested by IWiegeknann et al.l 
(l2006h . The vector components of the total magnetic force and 
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the total magnetic torque on the volume considered are given 
by six boundary integrals that must vanish if the magnetic field 
is force-free in t he full volume (Molodensky 1969; Aly 1984, 
fT989HLow|[T985h . The preprocessing changes the boundary val- 
ues of B within the error margins of the measurement in such a 
way that the moduli of the six boundary integrals are minimized. 
The resulting boundary values are expected to be more suit- 
able for an extrapolation into a force-free field than the original 
values. In the practical calculations, the convergence properties 
of the preprocessing iterations, as well as the calculated fields 
themselves, are very sensitive to small-scale noise and appar- 
ent discontinuities in the photospheric magnetograph data. This 
problem should, in principle, disappear if small spatial scales 
were sufficiently resolved. However, the numerical effort for that 
would be enormous. The small-scale fluctuations in the magne- 
tograms are also presumed to affect the solutio ns only in a very 
thin boundary layer close to the photosphere (Fu hrmann et alj 
12001) . Therefore, smoothing of the data is included in the pre- 
processing. 

The good performance of the optimization method, as in- 
dicated in Schrijver et al. (2006), encouraged us to develop 
a spherical versio n of the optimization code such as in 
IWiegelmannI (l2007h for a full sphere. In the first few sections of 
this paper, we describe a newly developed code that originates 
fr om a c artesian force-free optimization method implemented 
bv lWiegelm ann (2004). Our new code takes the curvature of the 
Sun's surface into account when modeling the coronal magnetic 
field in restricted area of the Sun. The optimization procedure 
considers six boundary faces, but in practice only the bottom 
boundary face is measured. On the other five faces, the assumed 
boundary data may have a strong influence on the solution. For 
this reason, it is desirable to move these faces as far away as 
possible from the region of interest. This, however, eventually 
requires that the surface curvature is taken into account. 

iDeRosa et al.l (l2009l) compared several nonlinear force-free 
codes in cartesian geo metry with stereos c opic reconstructed 
loops as produced by lAschwanden et al] ( l2008h . The codes 
used as input vector magnetograms from the Hinode-SOT-SP, 
which were unfortunately available for only a very small field 
of view (about 10 percent of the area spanned by STEREO- 
loops). Outside the Hinode FOV (field of view) line-of-sight 
magnetograms from SOHO/MDI were used and in the MDI- 
area, different assumptions about the transversal magnetic field 
have been made. Unfortunately, the comparison inferred that 
when different codes were implemented in the region outside 
the Hinode-FOV in different ways, the resulting coronal mag- 
netic field models produced by the separate codes were not con- 
sistent with the STEREO-loops. The recommendations of the 
authors are that one needs far larger high resolution vector mag- 
netograms, the codes need to account for uncertainties in the 
magnetograms, and one must have a clearer understanding of 
the photospheric-to-corona interface. Full disc vector magne- 
tograms will soon become available with SDO/HMl, but for a 
meaningful application we have to take the curvature of the Sun 
into account and carry out nonlinear force-free computations in 
spherical geometry. In this paper, we carry out the appropriate 
tests. We investigate first ideal model data and later data that con- 
tain artificial noise. To deal with noisy data and data with other 
uncertainties, we developed a preprocessing routine in spheri- 
cal geometry. While preprocessing does not model the details of 
the interface between the forced photosphere and the force-free 
base of the solar corona the procedure helps us to find suitable 
boundary conditions for a force-free modelling from measure- 
ments with inconsistencies. 



In this paper, we develop a spherical version of both the 
preprocessing and the optimization code for restricted part of 
the Sun. We follow the suggestion of Wiegelmann et al. (200^ 
to generalize their method of preprocessing photospheric vector 
magnetograms to spherical geometry just by considering the cur- 
vature of the Sun's surface for larger field of views. The paper 
is organized as follows: in Sect. 2, we describe an optimization 
procedure in spherical geometry; then, in Sect. 3, we apply it to a 
known nonlinear force-free test field and calculate some figures 
of merit for different boundary conditions. We derive force-free 
consistency criteria and describe the preprocessing procedure in 
spherical geometry in Sect. 4 and Sect. 5, respectively. In Sect. 
6, we use a known semi-analytic force-free model to check our 
method and in Sect. 5, we apply the method to different noise 
models. Finally, in Sect. 7, we draw conclusions and discuss our 
results. 

2. Optimization procedure 

Stationary states of the magnetic field configuration are de- 
scribed by the requirement that the Lorentz force be zero. 
Optimization procedure is one of several methods that have been 
developed over the past few decades to compute the most general 
class of those force-free fields. 

2. 7. Optimization principle in spliericai geometry 
Force-free magnetic fields must obey the equations 

(VxB)xB = 0, (1) 

V ■ B = (2) 

Equations ([T]) and (|2]i can be solved with the help of an op- 
timization principle, as proposed by IWheatland et alj ( l2000h 
and generalized by Wiegelmann ( 200-^ for cartesian geometry. 
The method minimizes a joint measure (La,) of the normalized 
Lorentz forces and the divergence of the field throughout the 
volume of interest, V. Here w e define a functional in spherical 
geometry ( IWiegelmannll2007h : 

Laj = J Cij(r,e,(l>)[B-^\('^ xB)xBf + \V ■ Bf'y^smedrdGdcp, 

(3) 

where (L>{r, 6, (p) is a weighting function and V is a computational 
box of wedge-shaped volume, which includes the inner physical 
domain V and the buffer zone(the region outside the physical 
domain), as shown in Fig.[3]of the bottom boundary on the pho- 
tosphere. The physical domain V is a wedge-shaped volume, 

o o 

with two latitudinal boundaries at ^min — 20 and flmax = 160 

o o 

, two longitudinal boundaries at (p^^m - 90 and 0n,ax = 270 , 
and two radial boundaries at the photosphere (r = 1^©) and 
r - 2Rq. The idea is to define an interior physical region V in 
which we wish to calculate the magnetic field so that it fulfills the 
force-free or MHS equations. We define V to be the inner region 
of V (including the photosphere) with o) - I everywhere includ- 
ing its six inner boundaries 6V' . We use the position-dependent 
weighting function to introduce a buffer boundary of nd - 6 
grid points towards the side and top boundaries of the compu- 
tational box, V. The weighting function, w is chosen to be con- 
stant within the inner physical domain V and declines to with 
a cosine profile in the buffer boundary region. The framed region 
in Figs.[3a-c) corresponds to the lower boundary of the physical 
domain V with a resolution of 48 x 62 pixels in the photosphere. 
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It is obvious that the force-free Eqs. ([T) and (|2| are fulfilled 
when La, equals zero. For fixed boundary conditions, the func- 
tional L„ in Eq. Q can be numerically minimized with the help 
of the iteration 

dB 

^^,F, (4) 

where // is a positive constant and the vector field F is calculated 
from 

F ^ojF + (Q„ X B) X Vw + (Qfe ■ B)Vw , (5) 

F = Vx(a„xB)-fiaX(VxB)+V(ni-B)-Q/,(V-B)+(Q;-HQ^)B, 

(6) 
(7) 

(8) 



n„ = B-^(V xB)xB, 
nb = B-^iV ■ B)B 



The field on the outer boundaries is always fixed here as 
Diric hlet boundary conditions. Relaxin g these boundaries is pos- 
sible dWiegelmann &Neukirchl [20031) and leads to additional 
terms. For w(r, 0, 0) = 1, the optimization method requires that 
the magnetic field is given on all six boundaries of V . This 
causes a serious limitation of the method because these data 
are only available for model configurations. For the reconstruc- 
tion of the coronal magnetic field, it is necessary to develop a 
method that reconstructs the ma gnetic field only from photo- 
spheric vector magnetograms (lWiegelmannll2004b . Since only 
the bottom boundary is measured, one has to make assumptions 
about the lateral and top boundaries, e.g., assume a potential 
field. This leads to inconsistent boundary conditions (see lAlvl 
Il989l regarding the compatibility of photospheric vector magne- 
tograph data). With the help of the weighting function, the five 
inconsistent boundaries are replaced by boundary layers and we 
consequently obtain more flexible boundaries around the phys- 
ical domain that will be adjusted automatically during the iter- 
ation. This diminishes the effect of the top and lateral bound- 
aries on the magnetic field solution inside the computational box. 
Additionally, the influence of the boundaries is diminished, the 
farther we move them away from the region of interest. 

The theoretical deviati on of the iterative Eq. (HI as outlined 
by IWheatland et al.l ( l2000l) does not depend on the use of a spe- 
cific coordinate system. Previous n umerical implementa tions of 
this method were demonstrated by IWiegeknannl ( |2007|) for the 
full sphere. Within this work, we use a spherical geometry, but 
for only a limited part of the sphere, e.g., large active regions, 
several (magnetically connected) active regions and full disc 
computations. Full disc vector magnetograms should become 
available soon from SDO/HMI. This kind of computational box 
will become necessary when the observed photospheric vector 
magnetogram becomes available for only parts of the photo- 
sphere. 

We use a spherical grid r, 0, (p with «,., «e, grid points in 
the direction of radius, latitude, and longitude, respectively. We 
normalize the magnetic field with the average radial magnetic 
field on the photosphere and the length scale with a solar radius. 
The method works as follows; 

1. We compute an initial source surface potential field in the 
computational domain from Br in the photosphere at r = 

IRQ. 

2. We replace Bg and B^ at the bottom photospheric boundary at 
r - IRq with the measured vector magnetogram. The outer 
radial and lateral boundaries are unchanged from the initial 
potential field model. For the purpose of code testing, we 
also tested different boundary conditions (see next section). 



3. We iterate for a force-free magnetic field in the computa- 
tional box by minimizing the functional L of Eq.© by ap- 

plying Eq.®. For each iteration step (k), the vector field F 
is calculated from the known field and a new field may 

simply be computed as = B**^' + f **^*Af for sufficiently 

small At. 

4. The continuous form of Eq.© ensures a monotonically de- 
creasing functional L. For finite time steps, this is also en- 
sured if the iteration time step dt is sufficiently small. If 
L{t + dt) > L{t), this step is rejected and we repeat this step 
with dt reduced by a factor of 2. 

5. After each successful iteration step, we increase dt by a fac- 
tor of 1.01 to ensure a time step as large as possible within 
the stability criteria. This ensures an iteration time step close 
to its optimum. 

6. The iteration stops if dt becomes too small. As a stopping 
criteria, we use dt < 10"^. 



2.2. Figures of merit 

To quantify the degree of agreement between vector fields B (for 
the input model field) and b (the NLFF model solutions) spec- 
ified on identical sets of grid points, we use five metrics that 
compare either local characteristics (e.g., vector magnitudes and 
directions at each point) or the global energy con tent in addition 
to the force and divergence integrals as defined inlS chrijver et al.l 
(|2006|) . The vector correlation (Cvec) metric generalizes the stan- 
dard correlation coefficient for scalar functions given by 



C,,c^J]Brb,/\J]\B,\^J]\b 



1/2 



(9) 



where B,- and 6, are the vectors at each point grid /. If the vector 
fields are identical, then Cuec = 1; if B, ± ft, , then C^ec - 0. 
The second metric, Ccs is based on the Cauchy-Schwarz 
inequalityda ■ b\ < \a\\b\ for any vector a and b) 

1 yr^ Bi b , 



where is the number of vectors in the field. This metric is 
mostly a measure of the angular differences between the vector 
fields: Ccs - 1, when B and b are parallel, and Ccs - ^1, if 
they are anti-parallel; Ccs - 0, if B, ± ft, at each point. 
We next introduce two measures of the vector errors, one nor- 
malized to the average vector norm, one averaging over relative 
differences. The normalized vector error E^j is defined as 



E^^Yj\^,-Bi\lY,\Bi\ 



The mean vector error /sm is defined as 



\Bi\ 



(11) 



(12) 



UnUke the first two metrics, perfect agreement between the two 
vector fields results inEM - Ei^ - 0. 

Since we are also interested in determining how well the models 
estimate the energy contained in the field, we use the total mag- 
netic energy in the model field normalized to the total magnetic 
energy in the input field as a global measure of the quality of the 
fit 

where e - 1 for closest agreement between the model field and 
the nonlinear force-free model solutions. 
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3. Test case and application to ideal boundary 
conditions 

3.1. Test case 

To te st the method, a kn own semi-analytic nonUnear solution is 
used. lLow & Lould 19901) presented a class of axisymmetric non- 
linear force-free fields with a multipolar character. The authors 
solved the Grad-Shafranov equation for axisymmetric force-free 
fields in spherical coordinates r, 0, and (p. The magnetic field can 
be written in the form 



B 



1 AdA^ dA^ V 



dr 



(14) 



where A is the flux function and Q represents the 0-component of 
B, depending only on A. The flux function A satisfies the Grad- 
Shafranov equation 



d^A l-ji-d-A ^dQ 



where ji - cos0. iLow & Loul (1 19901) derive solutions for 
dQ 

— = Q- = const, 
dA 

by looking for separable solutions of the form 



A(r, 0) 



(15) 



(16) 



(17) 



ILow & Loul (1 19901) suggested that these field solutions are 
the ideal solutions for testing methods of reconstruct- 
ing force-free fields from boundary values. They have 
become a standard test for nonlinear f orce-free e xtrap- 



olation codes in cartesian geometrv jAmari et al 



20061: IWheatland et al.l l2000t IWiegelmann & NeukirchI 



1999 



2003 



Yan & Lil 120061: llnhester & WiegelmannI 120061: ISchriiver etall 
20061) " because the symmetry in the solution is no longer obvi- 



ous after a translation that places the point source outside the 
computational domain and a rotation of the symmetry axis with 
respect to the domain edges. 

Here we use a Low and Lou solution in spherical coordi- 
nates. The solution used is labefled P i i with <1> = tt/IO, in Low 
& Lou's, notation (iLow &Loulll990h .' The original equilibrium 
is invariant in <p, but we can produce a (/(-variation in our coordi- 
nate system by placing the origin of the solution at / - 0.25 solar 
radius position from the sun centre. The corresponding configu- 
ration is then no longer symmetric in <f> with respect to the solar 
surface, as seen in the magnetic field map in the top row of Fig. 
[3] which shows the three components Br, Bg, and B^ on the pho- 
tosphere, respectively. We remark that we use the solution only 
for the purpose of testing our code and the equilibrium is not as- 
sumed to be a realistic model for the coronal magnetic field. We 
do the test runs on spherical grids (r, 0, (p) of 20 x 48 x 62 and 
40 X 96 X 124 grid points. 

3.2. Application to ideal boundary conditions 

Here we used different boundary conditions extracted from the 
Low and Lou model magnetic field. 

- Case 1: The boundary fields are specified on y'(all the six 
boundaries SV of V). 

- Case 2: The boundary fields are only specified on the photo- 
sphere (the lower boundary of the physical domain V). 



- Case 3: The boundary fields are only specified on the photo- 
sphere (the lower boundary of the physical domain V) and 
with boundary layers (at the buffer zone) of nd = 6 grid 
points toward top and lateral boundaries of the computa- 
tional box V. 

For the boundary conditions in case 1, the field line plot (as 
shown in Fig. [TJ agrees with original Low and Lou reference 
field because the optimization method requires all boundaries 
bounding the computational volume as boundary conditions. 

For the boundary conditions in case 2, we used an optimiza- 
tion code without a weighting function (nd - 0) and with a pho- 
tospheric boundary. Here the boundaries of the physical domain 
coincide with the computational boundaries. The lateral and top 
boundaries assume the value of the potential field during the it- 
eration. Some low-lying field lines are represented quite well 
(right-hand picture in Fig.[T]secondrow). The field lines close to 
the box center are of course close to the bottom boundary and far 
away from the other boundaries. The (observed) bottom bound- 
ary has a higher influence on the field here than the potential lat- 
eral and top boundary. Other field lines, especially high-reaching 
field lines, deviate from the analytic solution. 

For the boundary condition in case 3, we implemented an op- 
timization code with a weighting function of nd - 6 grid points 
outside the physical domain. This reduces the effect of top and 
lateral boundaries where B is unknown as a> drops from 1 to 
outward across the boundary layer around the physical domain. 

The comparison of the field lines of the Low & Lou model 
field with the reconstructed field of case 3 (the last picture in Fig. 
[T]i shows that the quality of the reconstruction improves signif- 
icantly with the use of the weighting function. Additionally, the 
size and shape of a boundary layer influences the quality of the 
reconstruction ( IWiegelmannll2004b for cartesian geometry. The 
larger computational box displaces the lateral and top bound- 
ary further away from the physical domain and its influence on 
the solution consequently decreases. As a result, the magnetic 
field in the physical domain is dominated by the vector magne- 
togram data, which is exactly what is required for application to 
measured vector magnetograms. A potential field reconstruction 
obviously does not agree with the reference field. In particular, 
we are in able to compute the magnetic energy content of the 
coronal magnetic field to be approximately correct. The figures 
of merit show that the potential field is far away from the true 
solutions and contains only 67.6% of the magnetic energy. 

The degree of convergence towards a force-free and 
divergence-free model solution can be quantified by the integral 
measures of the Lorentz force and divergence terms in the min- 
imization functional in Eq. ([3j, computed over the entire model 
volume V: 



X 



co{r, 0, (^)B"^|(V xB)xB\ sin 0drd0d(p, 

= I u){r,0,(l))\V ■ B^r^ s.vi\0drd0d^, 
Jv 



where Lj and Ld measure how well the force-free and 
divergence-free conditions are fulfilled, respectively. In Table 1, 
we list the figures of merit for our extrapolations results as intro- 
duced in previous section. Column 1 indicates the corresponding 
test case. Columns 2-4 show how well the force and solenoidal 
condition are fulfilled, where Col. 2 contains the value of the 
functional as defined in Eq.(3) and Lj and Ld in Cols. 3 and 4 
correspond to the first (force-free) and second (solenoidal free) 
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Table 1. Quality of our reconstructions with several figures of merit as explained in section 3.2. We compute the figures for the 
three different cases along with the model reference field and potential field. 



Model 




w 




II V - B lU 


WjxB lU 


Cvec 


Ccs 


En 


Em 


e 


Steps 


Time 


Original 


Spherical grid 20 x 48 x 62 
0.029 0.015 0.014 


1.180 


1.355 


1 


1 








1 






Potential 


0.020 


0.007 


0.014 


1.706 


1.091 


0.736 


0.688 


0.573 


0.535 


0.676 






Case 1 


0.006 


0.004 


0.002 


0.454 


0.774 


0.999 


0.983 


0.012 


0.016 


1.005 


10000 


7.14 min 


Case 2 


33.236 


7.806 


25.430 


47.843 


24.135 


0.757 


0.726 


0.397 


0.451 


0.745 


110 


1.28 min 


Case 3 


0.009 


0.006 


0.03 


0.367 


0.787 


0.994 


0.967 


0.187 


0.097 


0.989 


12011 


17.54 min 




Spherical grid 40 x 96 x 124 




















Original 


0.005 


0.003 


0.002 


0.38 


0.71 


1 


1 








1 






Potential 


0.30 


0.0003 


0.30 


0.44 


0.23 


0.67 


0.77 


0.70 


0.67 


0.75 






Case 1 


0.002 


0.001 


0.0006 


0.38 


0.32 


0.998 


0.999 


0.004 


0.007 


1.001 


12522 


Ih 21min 


Case 2 


26.27 


10.20 


16.07 


20.40 


30.53 


0.799 


0.759 


0.411 


0.456 


0.798 


5673 


Ih Imin 


Case 3 


0.24 


0.20 


0.04 


0.630 


0.747 


0.996 


0.971 


0.186 


0.112 


0.996 


12143 


4h 57min 



part of L„. The evolution of the functional \j x B\, and | V ■ B| 
during the optimization process is shown in Fig.|2] One can see 
from this figure that the calculation does not converge for case 2, 
because of the problematic boundaries where the fields are un- 
known. Column 5 contains the L°° norm of the divergence of the 
magnetic field 

||V-B|U=sup|V-B| 

xeV 

and Col. 6 lists the L°° norm of the Lorentz force of the magnetic 
field 

II J xB ||oo= sup|jxB|. 

xeV 

The next five columns of Table 1 contain diff'erent measurements 
comparing our reconstructed field with the semi-analytic refer- 
ence field. The two vector fields agree perfectly if Cvec, Ccs, and 
€ are unity and if and are zero. Column 12 contains the 
number of iteration steps until convergence, and Col. 13 shows 
the computing time on 1 processor 

A comparison of the original reference field (Fig.[Tla)) with 
our nonlinear force-free reconstructions (cases 1-3) shows that 
the magnetic field line plots agree with the original for case 1 
and case 3 within the plotting precision. Case 2 shows some de- 
viations from the original, but the reconstructed field lines are 
much closer to the reference field than the initial potential field. 
The visual inspection of Fig. [1] is supported by the quantitative 
criteria shown in Table 1 . For case 1 and case 3 the formal force- 
free criteria (L„, Lf, Lj) are smaller than the discretization error 
of the analytic solution and the comparison metrics show almost 
perfect agreement with the reference field. The comparison met- 
rics (of Table 1) show that there is discrepancy between the ref- 
erence field and case 2 as the magnetic field solution is affected 
by nearby problematic top and lateral boundaries. In Fig. [T] we 
compare magnetic field line plots of the original model field with 
a corresponding potential field and nonlinear force-free recon- 
structions with different boundary conditions (case 1 - case 3). 
The colour coding shows the radial magnetic field in the pho- 
tosphere, as also shown in the magnetogram in Fig. [3ja). The 
images show the results of the computation on the 20 x 48 x 62 
grid. 

4. Consistency criteria in spherical geometry 

A more fundamental requirement of the boundary data is its con- 
sistency with the force-free field approximation. As shown by 
[Molodensky ( 1969) and Aly ( 1989), a balance between the total 
momentum and angular momentum exerted onto the numerical 
box in cartesian geometry by the magnetic field leads to a set of 



boundary integral constraints on the magnetic field. These con- 
straints should also be satisfied on the solar surface for the field 
at the coronal base in the vicinity of a sufficiently isolated mag- 
netic region and in a situation where th ere is no rapid dynam ical 
development. As explained in detail in iMolodenskvl (|1974|) . the 
sense of these relations is that on average a force-free field can- 
not exert a net tangential force on the boundary or shear stresses 
along axes lying along the boundary. In summary, the boundary 
data for the force-free extrapolation should fulfill the following 
conditions: 

1 . The boundary data should coincide with the photospheric ob- 
servations within measurement errors. 

2. The boundary data should be consistent with the assumption 
of a force-free magnetic field above. 

3. For computational reasons (finite differences), the boundary 
data should be sufficiently smooth. 

Additional a-priori assumption is about the photospheric data are 
that the magnetic flux from the photosphere is sufficiently distant 
from the boundaries of the observational domain and that the net 
flux is balanced, i.e., 

J" B,.(r = IR„ 0, (l>)dD. = 0, (18) 

where S is the area of a bottom boundary of the physical domain 
on the photosphere. 

Generally, the flux balance criterion must be applied to the 
entire, closed surface of the numerical box. However, we can 
only measure the magnetic field vector on the bottom photo- 
spheric boundary and the contributions of the lateral and top 
boundary remain unspecified. However, if a major part of the 
known flux from the bottom boundary is uncompensated, the 
final force-free magnetic field solution will depend markedly 
on how the uncompensated flux is distributed over the other 
five boundaries. This would result in a major uncertainty on 
the final force free magnetic field configuration. We therefore 
demand that the flux balance is satisfied with the bottom data 
alone (|Wiegelmann & Inhester 2006). If this is not the case, we 
classify the reconstructio n problem a s not being uniquely solv- 
able within the given box. lAlvl (1 1 9891) used the virial theorem to 
define the conditions that a vector magnetogram must fulfill to 
be consistent with the assumption of a force-free field above in 
cartesian geometry. Here in this paper, we assume the force-free 
a nd torque-free c onditions for spherical geometry as formulated 
in lSakurail(ll994h . i.e.. 
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(a) Original 



(b) Potential 




(c) Case 1 



(d) Case 2 




(e) Case 3 

Fig. 1. The figure shows the original reference field, a potential field, and the results of a nonlinear force-free reconstruction with 
different boundary conditions (case 1-3, see text). The color coding shows B,- on the photosphere and the disc centre corresponds to 

o 

180 longitude. 



1 . The total force on the boundary vanishes 



^ [^(B^+B^-B;) sin Q sin 0-B,.Bfl cos Q sin 4>-BrB^ cos 0]t/Q = 0, 



(20) 



r [i(Bg+B^-B^)sin6lcos</i-B,.BgCos6icos0+B,.B0sin(^]<iQ = 0, ^ 

^ (19) J [-(B^ + B^-B;)cos0 + B,Besin0]jQ = O (21) 
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Fig. 2. Evolution of L„ (as defined in Eq. 3), max(\force\), and 
max{\dwB\) during the optimization process. The sohd line cor- 
responds to case 3, the dash-dotted line to case 1, the long- 
dashed line to case 2. 



2. The total torque on the boundary vanished 
J" B, {B^ cos 6 cos (f) + Eg sin ^)dQ. — 0, 

^ B, {B^ cos sin - Eg cos ^)dQ. — 0, 



(22) 

(23) 
(24) 



As with the flux balance, these criteria must in general, be ap- 
plied to the entire surface of the numerical box. Since we as- 
sumed that the photospheric flux is sufliciently concentrated in 
the center and the net flux is in balance, we can expect the mag- 
netic field on the lateral and top boundaries to remain weak and 
hence these surfaces do not represent a significant contribution 
to the integrals of the constraints above. We therefore impose the 
criteria on the bottom boundary alone. From this beginning, we 
use the following notation for simplicity: 



E-g^^-iBl + Bl-B'^), Es^ J^{Bj + Bl + Bl 



)dn. 



B\ - Bg COS 9 COS if> - B^ sin cf. 
Z?3 = Z?0 COS 6? COS <p + Bg sin ^ 



B2 - Bg COS sin + B^ cos , 
, B4 = B^ cos 9sin(p - Bg cos (p 



To quantify the quality of the vector magnetograms with respect 
to the above criteria, w e introduce three d i mensio nless param- 
eters similar to those in IWiegelmann et al ] (I2006h . but now for 
spherical geometry; 



1 . The flux balance parameter 



SfhLX 



See Appendix A for derivation of those torque-balance equations. 



Is \Br\dn 

2. The force balance parameter 

eforce-(j J" i^B sin 0COS - BrBi]c/Q| + | J" i^B sin 9 sin 
-BrB2]dQ.\ + \ J [Eg cos 6 + B,.Bg sin 6]dD.ljl Eg 

3. The torque balance parameter 

^torque 

=(| J B,.B3<iQ| + I J B,B4dQ.\ + \ J B,B^ sin GdQ.\)lEE 

An observed vector magnetogramis then flux-balanced and con- 
sistent with the force-free assumption if: <K 1, Sfon-e 1 
andeto;.^„e 1. 

5. Preprocessing method 

The strategy of preprocessing is to define a functional L of the 
boundary values of B, such that on minimizing L the total mag- 
netic force and the total magnetic torque on the considered vol- 
ume, as well as a quantity measuring the degree of small-scale 
noise in the boundary data, simultaneously become small. Each 
of the quantities to be made small is measured by an appropri- 
ately defined subfunctional included in L. The different subfunc- 
tionals are weighted to control their relative importance. Even 
if we choose a sufficiently flux balanced isolated active region 
(sfiux ^ 1), we find that the force-free conditions s force 1 and 
Etorque ^ 1 are not usually fulfilled for measured vector mag- 
netograms. We therefore conclude, that force-free extrapolation 
methods shoul d not be used directly on observed vector mag- 
netograms (see iGarvl (1200 lb for /? > 1 in photosphere), partic- 
ularly not on very noisy transverse photospheric magnetic field 
measurements. The large noise in the transverse components of 
the photospheric field vector, which is one order of magnitude 
higher than the LOS-field (~the transverse Bg and B^ at the bot- 
tom boundary), provides us freedom to adjust these data within 
the noise level. We use this freedom to drive the data towards 
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being more consistent with Aly's force-free and torque-free con- 
ditions. 

The preprocessing scheme of IWiegelmann eta n (I2006h m- 
volves minimizing a two-dimensional functional of quadratic 
form similar to the following: 

L - n\Lx + JU2L2 + A'3i3 + /^4i4 (25) 
Here we write the individual terms in spherical co-ordinates as: 

L\-iy ^ [Eg sin cos - B, B\\ sin + ^ ^ \E^^ sin sin 
p p 

-ByBx] sin of + V'b cos + BrBg sin 9] sin of , 

p 

(26) 

L2 = ( ^ BrB^ sin of + {^ BrB4 sin f)^ + ( ^ B,B^ sin^ of , 
p p p 

(27) 

L3 = {Br-Bi-obsf'+y^^ {Bg-Bgohsf + y^^ (B^-B^ohs)^, (28) 



L4 = ^ [(AB,)' + (ABe)' + (AB^f] (29) 



The surface integrals are replaced by a summation ^ — > 
lip sin dAOAcp, omitting the constant AdA(p over all p grid nodes 
of the bottom surface grid, with an elementary surface of 

sinOAcf) X AOy The differentiation in the smoothing term (L4) 
is achieved by the usual five-point stencil for the 2D-Laplace 
operator. Each of the constraints L„ is weighted by a yet unde- 
termined factor /Y„. The first term (« = 1) corresponds to the 
force-balance condition, and the next (« = 2) to the torque-free 
condition. The following term (« = 3) ensures that the optimized 
boundary condition agrees with the measured photospheric data, 
and that the last term (« = 4) controls the smoothing. The 2D- 
Laplace operator is designated by A. 

The aim of our preprocessing procedure is to minimize L so 
that all terms L„, if possible, become small simultaneously. This 
will yield a surface magnetic field: 



Bmin — argmin(L) 



(30) 



Besides a dependence on the observed magnetogram, the solu- 
tion in Eq.([25j now also depends on the coefficients fi„. These 
coefficients are a formal necessity because the terms L„ represent 
different quantities. By means of these coefficients, however, we 
can also give more or less weight to the individual terms in the 
case where a reduction in one term opposes a reduction in an- 
other This competition obviously exists between the observation 
term (« = 3) and the smoothing term (« = 4). The smoothing is 
performed consistently for all three magnetic field components. 

To obtain Eq.(l30]l by iteration, we need the derivative of L 
with respect to each of the three field components at every node 
(q) of the bottom boundary grid. We have, however, taken into 
account that Br is measured with much higher accuracy than Bg 
and B^. This is achieved by assuming that the vertical component 
is invariable compared to horizontal components in all terms 
where mixed products of the vertical and ho rizontal field com- 
ponents occur, e.g., within the constraints jWiegelmann et all 



I2OO6I) . The relevant functional derivatives of L are therefor^ 

dL 



d{Bg)„ 



dL 

d{B^\ 



-2fi 1 (Bg sin cos - B,. sin cos cos 0)^ x 

^ [E'g sin 6 cos (f) - B,B\] sin 6 
p 

+ 2// 1 (Bg sin^ sin - B, sin 6 cos 6 sin <f))q x 

^ [E'lj sin sin - B, B2] sin 6 

p 

+ 2nx{Bg sin 61 COS + B, sin- 6\x (31) 

^ [Eg COS + BrBg sin 0] sin 
p 

sin sin 0)^ ^ BrB3 sin 
/' 

- (Br sin cos (p)q ^ BrB4 sin flj 

p 

+ 2^3(Bg - Bgots)q + 2ll4(A(ABg)\, 

-2^1 (B^ sin^ cos (p + Br sin sin <f>)qX 

^ [E^j sin cos - B,.B 1 ] sin 

p 

+ 2/ii(B0sin^0sin0 - B,. sin0cos0)^x 

^ [Eg sin sin - BrB2] sin 

p 

+ 2yUi(B0 sin COS 0)q ^ [Eg cos + B,.Bg sin 0] sin 
p 

+ 2yU2[(Br COS COS (p sin 0)q ^ B,.B3 sin 

p 

+ (Br COS sin (f) sin 0)q ^ B,.B4 sin 

p 

+ (Br sin- 0)q BrB^ sin^ 0\ + 2fii(B^ - B^obs)q 
p 

+ 2^4(A(AB^))q, 
dL 



d(Br)q 



2fl3(Br - Br„bs)q + 2tl4(A(ABr))q 



(32) 
(33) 



The optimization is performed iteratively by a simple Newton or 
Landweber iteration, which replaces 



(Br)q^(B,)q-^i 
(Bg)q^(Bg)q-H 
{B^)q < (B^)q -II 



dL 
dL 



d(Bg)q' 

dL 
d(B^)q 



(34) 
(35) 
(36) 



at every step. The convergence of this scheme towards a solution 
of Eq.(|25]l is obvious: L has to decrease monotonically at every 
step as long as Eqs.(l26]l-(l28ll have a nonzero component. These 



^ See Appendix B for partial derivative of L4 with respect to each of 
the three field components. 
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Fig. 3. Top row: vector magnetogram derive from the Low and Lou solution. From left to right the three components B, , Bg & B^ are 
shown). Middle row: the same magnetogram as in the first row, but with noise added (noise model I). Bottom row: magnetogram 
resulting from preprocessing of the disturbed magnetogram shown in the second row. The magnetic fields are measured in gauss. 
The vertical and horizontal axes show latitude, 6 and longitude, on the photosphere respectively. 



terms, however, vanish only if an extremum of L is reached. 
Since L is fourth order in B, this may not necessarily be a global 
minimum; in rare cases, if the step size is handled carelessly, 
it may even be a local maximum. In practical calculation, this 
should not, however, be a problem and from our experience we 
rapidly o btain a minimum B,ni„ of L, once the parameters //„ are 
specified dWiegeknann et al]|2006l) . 

6. Application to different noise-models 

We extract the bottom boundary of the Low and Lou equilibrium 
and u se it as input for our extrapolation code (see Wiegelmann 
|2004|) . This artificial vector magnetogram (see first row of Fig. [3) 
extrapolated from a semi-analytical solution is of course in per- 
fect agreement with the assumption of a force-free field above 
(Aly-criteria) and the result of our extrapolation code was in rea- 
sonable agreement with the original. True measured vector mag- 
netograms are not ideal (and smooth) of course, and we simulate 
this effect by adding noise to the Low and Lou magnetogram 



dWiegelmann et al]|2006l) . We add noise to this ideal solution in 
the form: 

Noise model I: 

6Bi - rii ■ r„ ■ y/B], where «/ is the noise level and r„ a random 
number in the range -L...1. The noise level was chosen to be 
«/ = 10.0 for the transverse magnetic field (Bg, B^) and «/ = 0.5 
for B, . This mimics a real magnetogram (see the middle row of 
Fig.[3j with Gaussian noise and significantly higher noise in the 
transverse components of the magnetic field. 

Noise model II: 
6Bi = III ■ r„, where «/ is the noise level and r„ a random number 
in the range The noise level was chosen to be «/ = 20.0 

for the transverse magnetic field {Bg, B^) and «/ =1.0 for B, ). 
This noise model adds noise, independent of the local magnetic 
field strength. 

Noise model III: 
5B,- - constant, 6B, - f^'""" , where we choose a constant 

noise level 6B, of 1 and a minimum detection level 6Bt„,m - 20. 
This noise model mimics the effect in which the transverse 
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(c) Field from noisy data (d) Field from preprocessed data 



Fig. 4. a) Some field lines for the original Low and Lou solution, b) Potential field reconstruction, c) Nonlinear force-free reconstruc- 
tion from noisy data (noise model 1) without preprocessing, d) Nonlinear force-free reconstruction from noisy data (noise model I) 
after preprocessing the vector magnetogram with our newly developed spherical code. 



Table 2. Figures of merit for the three different noise models with and without preprocessing along with model reference field and 
potential field. 



Model 


Preprocessed 


L 


U 


Li 


IIV-BIU 


II jxB lU 


c 


Ccs 


En 


Em 


e 


Steps 


Original 




0.029 


0.015 


0.014 


1.180 


1.355 


1 


1 








1 




Potential 




0.020 


0.007 


0.014 


1.706 


1.091 


0.736 


0.688 


0.573 


0.535 


0.676 




Noise model 1 


No 


22.015 


8.612 


13.403 


25.531 


11.671 


0.819 


0.767 


0.337 


0.421 


0.861 


1337 


Noise model 1 


Yes 


0.105 


0.066 


0.039 


1.746 


1.806 


0.951 


0.947 


0.197 


0.105 


0.964 


12191 


Noise model 11 


No 


18.957 


7.915 


11.042 


23.089 


9.871 


0.828 


0.774 


0.321 


0.417 


0.869 


1484 


Noise model 11 


Yes 


0.097 


0.057 


0.040 


1.533 


1.617 


0.963 


0.951 


0.191 


0.099 


0.971 


11423 


Noise model 111 


No 


17.718 


7.615 


10.103 


20.763 


8.992 


0.859 


0.781 


0.310 


0.402 


0.873 


1497 


Noise model 111 


Yes 


0.081 


0.043 


0.038 


1.382 


1.407 


0.979 


0.957 


0.189 


0.098 


0.982 


10378 



noise level is higher in r egions of low magnetic field strength 
dWiegelmann et al.ll2006l) . 

The bottom row of Fig. [3] shows the preprocessed vector 
magnetogram (for noise model I) after applying our procedure. 
The aim of the preprocessing is to use the resulting magnetogram 
as input for a nonlinear force-free magnetic field extrapolation. 

Figure|4]shows in panel a) the original Low and Lou solution 
and in panel b) a corresponding potential field reconstruction. In 
Fig. |4] we present only the inner region of the whole magne- 
togram (marked with black rectangular box in Figl3a)) because 
the surrounding magnetogram is used as a boundary layer (6 grid 



points) for our nonlinear force-free code. The computation was 
done on a 26 X 60 X 74 grid including a 6 pixel boundary layer to- 
wards the lateral and top boundary of the computational box V. 
In the remaining panels of Fig. H] we demonstrate the effect of 
the noise model (I) on the reconstruction. The noise levels were 
chosen so that the mean noise was similar for all three noise 
models. Fig. |4] (c) shows a nonlinear force-free reconstruction 
with noisy data (noise model I, magnetogram shown in the cen- 
tral panel of Fig.[3]l, and Fig.|4](d) presents a nonlinear force-free 
reconstruction after preprocessing (magnetogram shown in the 
bottom panel of Fig. [3]l. After preprocessing(see Fig. |4]d), we 
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Fig. 5. Vector correlation plotted against noise level for noise 
model I. 



achieve far closer agreement with the original solution (Fig. |4] 
a). Field lines are plotted from the same photospheric footpoints 
in the positive polarity. 

For the other noise models II and III, we find that the pre- 
processed data agree more closely with the original Fig.|4](a). 
We check the correlation of the original solution with our recon- 
struction with help of the vector correlation function as defined 
in®. 

Table 2 confirms the visual inspection of Fig. |4] The cor- 
relation of the reconstructed magnetic field with the origi- 
nal improves significantly after preprocessing of the data for 
all noise models. We knew already from previo us studies 
jWiegelmann & NeukirchI 120031: IWiegelmannll2004l) that noise 
and inconsistencies in vector magnetograms have a negative in- 
fluence on the nonlinear force-free reconstruction, and the pre- 
processing routine described in this paper shows us how to over- 
come these difficulties in the case of spherical geometry. As 
indicated by Fig. [5] the higher the noise level we have added 
to the original magnetogram, the smaller the vector correlation 
will be for the field reconstructed from the magnetogram with 
noise, compared with the reference field. However, the corre- 
sponding vector correlations for the field reconstructed from the 
preprocessed magnetogram has no significant change as the code 
largely removes the noise we have added to the original magne- 
togram with diff'erent noise levels. 



boundary layer and reaches zero at the boundary of the com- 
putational box. At the boundary of the computational box, we 
set the field to have the value of the potential field computed 
from Br at the bottom boundary. Our test calculations show that 
a finite-sized weighted boundary yields far more reliable results. 
The depth nd of this buffer boundary influences the quality of 
reconstruction, since the magnetic flux in these test cases is not 
concentrated well inside the interior of the box. 

In this work, we have presented a method for preprocess- 
ing vector magnetogram data to be able to use the preprocess- 
ing result as input for a nonlinear force-free magnetic field ex- 
trapolation with help of an optimization code in spherical ge- 
ometry. We extended the preprocessing routine developed by 
IWiegelmann et al.l (|2006|) to spherical geometry. As a first test 
of the method, we use the Low and Lou solution with added 
noise from different noise models. A direct use of the noisy pho- 
tospheric data for a nonlinear force-free extrapolation showed no 
good agreement with the original Low and Lou solution, but af- 
ter applying our newly developed preprocessing method we ob- 
tained a reasonable agreement with the original. The preprocess- 
ing method changes the boundary data within their noise limits 
to drive the magnetogram towards boundary conditions that are 
consistent with the assumption of a force-free field above. The 
transverse field components with higher noise level are modified 
more than the radial components. 

To carry out the preprocessing, we use a minimization prin- 
ciple. On the one hand, we control the final boundary data to 
be as close as possible (within the noise level) to the original 
measured data, and the data are forced to fulfill the consistency 
criteria and be sufficiently smooth. Smoothness of the boundary 
data is required by the nonlinear force-free extrapolation code, 
but also necessary physically because the magnetic field at the 
basis of the corona should be smoother than in the photosphere, 
where it is measured. In addition to these, we found that adding a 
larger amount of noise to the magnetogram decreases its vector 
correlation with the model reference field whenever we recon- 
struct it without preprocessing. 

We plan to use this newly developed code for future mis- 
sions such as SDO (Solar Dynamic Observatory) when full disc 
magnetogram data become available. 
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7. Conclusion and outlook 

In this paper, we have developed and tested the optimization 
method for the reconstruction of nonlinear force-free coronal 
magnetic fields in spherical geometry by r estricting the code t o 
limited parts of the Sun, as suggested by IWiegelmannI (|2007|) . 
The optimization method minimizes a functional consisting of 
a quadratic form of the force balance and the solenoidal con- 
dition. Without a weighting function, all the six boundaries are 
equal likely to influence the solution. The effect of top and lat- 
eral boundaries can be reduce d by introducing a boundary layer 
around the physical domain d WiegelmannI 120041) . The physical 
domain is a wedge-shaped area within which we reconstruct the 
coronal magnetic field that is consistent with the photospheric 
vector magnetogram data. The boundary layer replaces the hard 
lateral and top boundary used previously. In the physical domain, 
the weighting function is unity. It drops monotonically in the 
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Appendix A: Torque-balance equations 

To solve the torque-balance equations in Eqs.(l22b-(|24|). we as- 
sume that the volume integral of torque in the computational box 
vanishes to fulfill the force-free criteria. 



X 



(r X F)dV = 0, 



where the force F is 



f = — (VxB)xB 



(A.l) 



(A.2) 



Substituting Eq.( IA.2l i into Eq. dA.ll ) and using vector identity 

V(A-B) = Ax(VxB)+Bx(VxA) + (A-V)B + (B- V)A, (A.3) 

along with the Gauss divergence theorem we can find the fol- 
lowing expression 

1 (rxf)t/y=i 1 B^rxdS)- f (rxB)(B £/S) = 0. (A.4) 
Jv 2 

The vector dS is directed into the volume V, with dS = 



R%. sin 6d0d<per- The origin of the vector r = RQCr is taken to be 



the centre of the Sun. Therefore, Eq.( |A.4t reduces to the form 
(r X B){B ■ dS) = 0, (A.5) 



I' 



since (r x dS) - 0. Where B - B,e, + BgCg + B^e^, we have 
B dS ^ R^B,. sin 6d6d(f) - R^BrdQ.. By substituting the vector 
rxB into Eq. (IA.5b . one can find 



^ [ - RQB^ee + RoBge^jB^Ry^l = 0. (A.6) 

Hence the torque balance along the A:-component will be 

e.v [ ^ ( - ReB^ee + RQBge^)B,.Ryn] = 0, (A.7) 

where e ^ is the unit vector along the x-axis. Normalizing this 
equation by assuming that Rq = 1, Eq. (IA.7b reduces to 



X 



BriB^cosecoscp + Bg sin (l))dQ. = 0, (A.8) 



Similarly one can obtain Eqs. ( |23] | and ( l24l i. 



Appendix B: Partial derivative of L4 

We derive the partial derivative of L4 with respect to each of 
the three magnetic field components in its discretized form as 
indicated in Eqs.(|3TTi-(|33Tl. We used a five-point stencil on the 
photospheric boundary for Laplace in L4. Those derivatives are 
carried out at every node (q) of the bottom boundary grid. The 
partial derivative of Eq. ( |29] l with respect to B,-, for instance can 
be written as 

We demonstrated the effect of the derivative by using the conven- 
tional Laplacian ABr in one dimension using three-point stencil 
with geometry-dependent coefficients c & a. Then 

(AB,.)„ = a{Br)p-i + c{B,)p + a(Br)p^x , (B.2) 
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and after substituting Eq.( |B.2t into the derivative term in 
Eq. dB.lt . we find 



d(B,), diB,),^"'""" " (B.3) 

-a6p-\^q + C6p,q + aSp+l^q 

Therefore, using equation Eq. dB.Sb . we can reduce Eq. dB.lb to 
dU 



=2 ^{lS.Br)p{a6p-\_^q + C5p,q + a6p+l,q) 

p 

=2 ^ [fl(AB,)^+i + c(AB,), + a(AB,)^-i] (b.4) 
p 

=2;^(A(A5,))^. 



One can similarly derive the partial derivative of L4 with respect 
to the other two field components. 



